% To do likelihood ratio test
clear all;
clc;
disp('To do likelihood ratio test_MK1');
load Counts_TP06YRJfinR5G2501000kMR7.mat;
% load Avg_Prob_TP06YRJfinR5G1000kMR7.mat;
P1=zeros(7,7);

% take the needed data from count. totaly 12yrs from 1998-2009
N1(:,:)=count(:,3:9,3:9);
N1(:,8)=count(:,3:9,11);
ttldata=sum(N1(:));
filename='TP06_Y_RJ_fin_R5_G250_1000k_MR7_Likelihood Ratio Test_ND_MK1.xls';

% to calculate log (LND) hypothesis 
N0=sum(N1(1:7,:));
P0=N0(:,:)./sum(N0(:,1:8));
%XX(:,:)=log(P0(:,:)).*N0(:,:)
XX(:,:)=log(P0(:,:)).*N0(:,:);
LLND=sum(XX(:));

% to calculate log (LMK1) hyphothesis
for i=1:7
    for j=1:8
        P1(i,j)=N1(i,j)./sum(N1(i,1:8));
    end
end
%X(:,:)=log(P1(:,:)).* N1(:,:)
X(:,:)=log(P1(:,:)).* N1(:,:);
LLMK1=sum(X(:));

%Likelihood Ratio Test
invalid_Num_N0 = prod(size(find(N0==0)));
invalid_Num_N1 = prod(size(find(N1==0)));
dof = prod(size(P1)) - prod(size(P0)) - invalid_Num_N0 - invalid_Num_N1;% minus the number of zero element in P0 and P1.
[LRh,LRp,LRstat,cV] = lratiotest(LLMK1,LLND,dof);
 
sheetlist1=strcat('Prob of j');
sheetlist2=strcat('Prob of ij');
sheetlist3=strcat('ObservationNum N0');
sheetlist4=strcat('ObservationNum N1');
sheetlist5=strcat('Parameters');
xlswrite(filename, P0, sheetlist1);
xlswrite(filename, P1, sheetlist2);
xlswrite(filename, N0, sheetlist3);
xlswrite(filename, N1, sheetlist4);
d = {'LRh','LRp','LRstat','cV','dof','LLND','LLMK1','ttlData'; LRh, LRp, LRstat,cV,dof,LLND,LLMK1,ttldata};
xlswrite(filename, d, sheetlist5)


